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Abstract 

Individual locations of many neuronal cell bodies (> 10^) are 
needed to enable statistically significant measurements of spatial or- 
ganization within the brain such as nearest-neighbor and microcolum- 
narity measurements. In this paper, we introduce an Automated Neu- 
ron Recognition Algorithm (ANRA) which obtains the (x,y) location 
of individual neurons within digitized images of Nissl-stained, 30 mi- 
cron thick, frozen sections of the cerebral cortex of the Rhesus mon- 
key. Identification of neurons within such Nissl-stained sections is 
inherently difficult due to the variability in neuron staining, the over- 
lap of neurons, the presence of partial or damaged neurons at tissue 
surfaces, and the presence of non-neuron objects, such as glial cells, 
blood vessels, and random artifacts. To overcome these challenges 
and identify neurons, ANRA applies a combination of image segmen- 
tation and machine learning. The steps involve active contour seg- 
mentation to find outlines of potential neuron cell bodies followed 



1 



by artificial neural network training using the segmentation prop- 
erties (size, optical density, gyration, etc.) to distinguish between 
neuron and non-neuron segmentations. ANRA positively identifies 
86 ± 5% neurons with 15 ± 8% error (mean ± st.dev.) on a wide 
range of Nissl-stained images, whereas semi-automatic methods ob- 
tain 80 ± 7%/17 zb 12%. A further advantage of ANRA is that it 
affords an unlimited increase in speed from semi-automatic methods, 
and is computationally efficient, with the ability to recognize ^100 
neurons per minute using a standard personal computer. ANRA is 
amenable to analysis of huge photo- montages of Nissl-stained tissue, 
thereby opening the door to fast, efficient and quantitative analysis 
of vast stores of archival material that exist in laboratories and re- 
search collections around the world. The ANRA software is available 
at http : //physics . bu . edu/^^ainglis/ANRA/ , 



1 Introduction 



Since the 1980s, the application of unbiased stereological approaches to quan- 
tify objects of biological interest has allowed for rigorous measurements of 
many parameters of brain structure including total neuron number, area, and 
volume. These approaches are based on systematic random sampling from 



defined regions of interest using unbiased estimators (May hew, 1991; Schmitz 



& Hof, 2005). While these measurements have produced extremely valuable 



insights into the structural organization of the brain, including age-related 



preservation of neuron numbers (Peters et a/., 1998), these "first order" stere 



ological parameters only partially describe the structural organization of the 
brain, as they cannot efficiently quantify "second order" parameters that 
measure more complex spatial properties of neuron organization, such as the 



nearest neighbor arrangement ( 


Asare 


, 19961 


Schmitz et al , 


20021 Duyckaerts 


& Godefroy, |2000l |Krasnoperov & Stoyan, 2004t Hof et al\ |2003l |Urbanc 


et al , 


2002|) and arrangement into mini- or microcolumns ( 


Cruz et a/., 2005; 



Buldyrev et a/., 2000; Buxhoeveden & Lefkowitz, 1996). 



Several approaches can be used to quantify "second order" parameters. 



Stereological methods can quantify nearest-neighbor arrangement (Schmitz 



et al, 2002), but the methods are labor intensive and would be difficult to 
apply to large brain areas. Image Fourier methods do not require manual 
marking of neuron locations and can quantify "vertical bias" of objects within 



an image ( jCasanova et al\ [2006|) , but do not discern between the contribution 
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from glial and neuronal cell bodies. 

Alternatively, pair correlation methods use concepts from statistical physics 



to calculate correlation properties such as ID nearest-neighbor (Urbane et al 



2002) and 2D microcolumnar organization (Buldyrev et a/., 2000) of neurons, 



as well as more discerning properties of spatial arrangement, such as the 



strength of microcolumnar order and microcolumnar width and length (Cruz 



et a/., 2005). The multitude of spatial organization quantities that can be 



calculated with pair correlation analysis makes it appealing to apply to large 
brain areas. To do that, we first need to address the major challenge to this 
approach: how to obtain the necessarily large number of neuron locations 
(10^ — 10^ locations per measurement) to get statistically significant results 
(see Sec. 2.7 and Discussion) over large regions of the brain, reaching ^ 10^ 
for a large study. The acquisition of such numbers of neurons by manually 
or semi-automatically identifying and marking the location of each is pro- 
hibitively time-consuming and open to user bias. Hence, correlative analysis 
of spatial relationships among neurons (as well as non-stereology based cell 
counts (Todtenkopf et a/., 2005] )) would be dramatically facilitated by an 
automatic method for identifying and locating the visible centers of neurons 
accurately and efficiently. 

While various other immunhistochemical methods could facilitate auto- 
mated discrimination of neurons and glia better than Nissl, there are im- 
portant advantages to develop automated methods for Nissl-stained tissue. 
Nissl-staining is the least expensive, easiest applied method for staining both 
neurons and glia. Furthermore, there are thousands of unique and often irre- 
producible collections of Nissl-stained brain material in clinical and research 
labs around the world that could be analyzed using the ANRA. 

There are several challenges to automatically retrieve neuron locations 
from two-dimensional digitized images of Nissl-stained brain tissue (Fig. [2^). 
A major challenge is to distinguish between neuron and non-neuron objects, 
including staining errors, tissue folds, and dirt particles, as well as blood 
vessels and glial cells. Another challenge is to identify neurons that differ 
almost as widely from each other as they do from non-neuronal objects. 
Neuron cell bodies are naturally diverse in size and shape and have different 
orientations with respect to their dendrite and axon processes. Neurons can 
also be cleaved at the cutting surface or damaged by the cutting process, 
which affects their shape in the tissue. These variables lead to diverse neuron 
cell profiles within the tissue slice. A further challenge is to discriminate 
between neurons that overlap, a common finding as tissue sections are 3D 
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volumes projected onto a 2D image. 

There are currently several published approaches to automatic retrieval 
of cell bodies from images. Some methods use segmentation techniques based 



on thresholding (Slater et al, 1996 



et al., 2003), watershed (Lin et al. 



Benah et al, 2003), Potts model (Peng 



2005), and active contours (Ray et al. 



2002 ) . Others use trained neural networks to mark appropriately sized "pixel 



patches" as cells of interest. The "pixel patch" training methods use artificial 



neural networks (Sjostrom et a/., 1999), local linear mapping (Nattkemper 



et a/., 2001), Fischer's linear discriminant (Long et al., 2005), and support 



vector machines (Long et al., 2006). Another method based on template 



matching has been recently introduced by Costa & BoUt (2006). 

In this paper we introduce and test an Automatic Neuron Recognition 
Algorithm (ANRA) (Fig. [T]) which uses a combination of segmenting and 
training to overcome the challenges of retrieving neuron location in Nissl- 
stained tissue sections. ANRA automatically identifies neurons from digital 
images and retrieves their (x,y) locations. 



2 Methods 

2.1 Image Input and Preprocessing 

The inputs for ANRA are photomicrographs of 30 micron thick Nissl-stained 
tissue section taken at lOx magnification and a resolution of 1.5 microns per 
pixel. Because the 30 micron tissue section shrinks during processing to a 
thickness of less than 10 microns, all of the tissue is in focus when viewed 
at microscopic magnifications of 20X or lower, thus the 2D image properly 
represents neuron locations . Since the color information is not as useful in 
the monotone Nissl-stained images (Fig.|2^) the images are converted to gray 
scale images ranging from (black) to 255 (white). 

The photomicrographs are taken from different areas of the brain from dif- 
ferent subjects at different times. Therefore, images are of different "quality" , 
refiecting a combination of variations in morphology, staining, slide prepa- 
ration, and digitization (Fig. ^p)- To reduce this variability, the images are 
first "normalized" such that every image has the same background and fore- 
ground average optical density. This is done by thresholding each image into 
foreground and background pixels and finding the average optical density for 
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the foreground and background separately. For each image, the optical den- 
sity histogram is then shifted to match the foreground/background averages 
of an ideal image (Fig. [3]a). Fig. shows the images final normalization as 
compared to the original images in Fig. [2|3. This preprocessing step removes 
most of the image variations due to processing (staining, slide preparation, 
digitization, etc.) and is a key step toward applying ANRA to an unlim- 
ited number of images that do not vary drastically in intrinsic morphological 
differences (neuron density, shape, size, etc.). There is no need for other 
preprocessing steps such as blurring or sharpening since ANRA, by design, 
overcomes the challenges of noisy images and weak boundary information. 



2.2 Main segmentation tool: OSM 

Here we describe the segmentation procedure presented in Fig. [T} called the 
overall segmentation method (OSM). 



2.2.1 Over-marking the image 

An initial step of the segmentation process is "seeding" the image with one 
or more points for each possible neuron cell body. A combination of two 
methods is used (Fig. [i^): a hexagonal grid of points is placed over the 
thresholded foreground of the image and the center points of objects identified 
by the traditional watershed segmentation (Javi, 2002). 



2.2.2 Active contour segmentation 



We employ active contour segmentation with statistical shape knowledge (Cre- 



mers et a/., 2000) because the method is designed to overcome the challenges 



of noisy images and missing boundary data, the main identification chal- 
lenge in Nissl-stained tissue. Also, the method uses low-dimensional shape 
representations which are ideal for modeling cell contours (outlines of cells). 
Because the image is initially over-marked, the calculations of contour split- 
ting (Zimmer et a/., 2002) are not needed. 



The image fij is a digital image of sliced brain tissue which defines the 
optical density (gray scale value) of each pixel (ij). We assume that the 
image contains at least one type of object of interest (neurons) mixed with 
other objects (non-neurons). The goal of a single run of the segmentation 
is to "segment" a single object of interest (a single neuron) from the rest of 
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the image (all other neurons, non-neurons, and background). It does this 
by "evolving" a loop of pixels called a contour (C) from a circle of typical 
neuron diameter (12/xm) starting at one of the over-marked starting points, 
to a location and shape that surrounds a potential neuron ceU body (Fig.|4|b). 
This process is repeated for each starting location until all starting locations 
have been exhausted. 

The movement of C is controlled by a set of points called control points 



{(^715 yn)}n=i N which wc usc the compact notation (Cremers et a/., 2000) 



z = (ri, Tat) = (xi, yi, Xat, Vn) (1) 

The control points are parameters in a closed quadratic Bezier-spline (B- 
spline) curve (Blake fc Isard[ |1998 ) that define the exact location (pixels) 



of C (see Fig. [5] for definition). Hence, C moves and changes shape by the 
iterative motion of the control points z. At each time step, each control point 
z makes a small movement towards encircling an object close to its starting 
location by minimizing a total energy E based on two energy considerations. 
Ems and Ec'. 



E{f,u,C) = EMs{f.u) + aE,{C) . (2) 

A qualitative understanding of the energy terms is presented in Fig. |6| Ems 
is the Mumford-Shah energy term, which determines how well the contour 
separates lighter and darker gray scale regions in the image fij. Ec{C) is 
the contour energy term, which quantifies the similarity of the contour to a 
previously chosen set of training shapes (in our case, the training shapes are 
oval-like). Ems is high when C does not separate different contrasts well, 
and is low if it does. Ec{C) is high if the shape is very contorted, and low 
if it is oval-like, a changes the relative infiuence of the two energy terms. If 
a is a high value, then C will evolve into a rigid perfect oval, ignoring all 
image information. If a is zero, then C will surround any nearby object in 
the image with no regard to the final shape of C . When the two energy terms 
are balanced with an appropriate a and the system is evolved to minimize 
E then objects in an image are encircled properly. Fig. [7] shows a typical 
evolution of C with an appropriate a value. IS a variable image, similar to 
a blurred version of /^j, which is used in the algorithm, as described below. 

The Mumford-Shah energy term EMs{f) '^) quantifies the alignment of the 
contour with edges in the image fij\ 
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EMs{f,u) = -Y.{{hj-un{t)r + X'\Vun{t)\'} (3) 

where A is the Mumford-Shah energy parameter that determines relative 
strengths of the terms. \Vuij{t)\'^ is the square of the magnitude of the 
picture gradient: 



\Vuij{t)\' 



du 
dx 



' fduV 



[Ui+l,j{t) - Ui_ij{t)f + Kj+l(t) - Ui^j-l{t)Y 



(4) 

It should be noted that Cremers et al. ( 2000D includes an additional term 
i^\\C\\ to Eq. |3| which minimizes the length ||C|| of the contour within its 
evolution. We do not include this term because it adds an additional free 
parameter and does not contribute to the functionality of the algorithm when 
identifying cell shaped objects. 

Eq. |3] is differentiated with respect to control point movement. Setting 
the solution of the differentiation to a minimum of EMs{f^u{t)) gives the 
evolution equation for each individual control point n = 1..N during each 



iteration dt (Mumford fc Shah, 1989): 



in (t) = (e+ - 
yn{t) = (e+ -e")n. 



(5) 



where e+ and e~ are Ems (Eq. [3]) summed over the single line of pixels right 
outside (e+) and right inside (e~) the segment of C centered around control 
point (xn^Vn) (Fig. [s]). n^; and n^; are the outer normal vectors of C at 
each control point in the x and y direction respectively, x = dx/dt and 
y = dy/dt ^ where t is the artificial time parameter. 

Eq. [3]is then differentiated with respect to the variable image Uij. Setting 
the solution to a minimum of EMs{f^u{t)) gives the evolution equation for 



each pixel Uij during each iteration dt (Mumford & Shah, 1989): 



Uij {t + dt) 



^^J(^) + {k - ^u-(i) + y^^u,j{t)} dt if ij 3 C 



U,j{t) 



if ij G C 



(6) 



At t = 0, tijj(O) = fij. V'^Uijit) is the Laplacian in 2-D Cartesian coordinates: 
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Eq.|6| describes a diffusion {V Uij{t)) process limited by the original image 
[fij — Uij{t)). The key component is that Uij never evolves at the pixels 
that make up C . Uij becomes stable once C separates contrasted regions. 
Therefore, minimizing Ems tends to evolve C so that the gray scale values 
vary slowly (smoothly) in the areas inside and outside the contour but vary 
strongly (discontinuously) across the contour C. 

The contour energy term Ec affects the shape of the contour irrespective 
of the images fij and Uij. Ec is minimized for contour shapes most similar 
to a previously chosen set of training shapes x = {^i^ ^2, •••}• The energy is 
calculated using the following equation: 

E,{C) = ]^{z-zof E-^{z-zo) , (8) 

where the vector and the matrix E (with an inverse contain the mean 
and covariant information of the previously chosen set of training shapes 
X = {^1,^2,...}: 

^0 = (9) 

E = {{z,-z^f{z,-Zo)) , (10) 

Here <> denotes the sample average. is a 2N vector and i7 is a 2Nx2N 
matrix. Creating zq and U for a set of shapes x = {^i, ^2, •••} is equivalent to 



modeling the distribution of shapes in M as a Gaussian distribution (Cre- 



mers et a/.[|2000[ ). 

To minimize Ec{C)^ the following evolution equation for each control 
point is used: 

z{t) = {z{t) - zo) . (11) 

Combining the two equations [5] and [TT] gives the final evolution equation 
for each control point n during each iteration: 



Xn {t + dt) = Xn (t) + {(e+ - e") + a [IJ-^ {z{t) - ^o)]2n-i} dt 
Vn {t + dt) = Vn {t) + {(e+ - e") riy + a [S-^ {z{t) - zo)]^J dt , 



(12) 
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recalling that e+ and e~ are Ems (Eq. [s]) summed over the single line 
of pixels right outside (e+) and right inside (e~) the segment of C centered 
around control point (x^, yn) (Fig. Isl), and are dependent on A and Uij. 



The evolution of the contour is driven by Eqs. (12 6), with variables Uij 
and contour points (xi, ^i, Xa^, Vn)- Note that Eqs. ( T2[6 ) are coupled and 
must be solved simultaneously. 



Uij 



Performing a step by step evolution of the control points (Eq. [12]) and 
(Eq. |6]), C evolves in the following way: If C begins to change into a 
contorted, non-ovular shape to minimize Ems (such as "leaking" out of an 
area of weak or missing boundary information in the image), then E^ will 
increase, hence there will be a force opposing the movement. Similarly, if 
the contour begins to move back to a perfect oval to minimize £'c. Ems will 
increase and thus limit such a change. When a local minimum is reached and 
the contour no longer moves, the points internal to the contour are saved, 
and the process starts again at a new location until all starting locations are 
exhausted. 

There are several free parameters (a. A, rft, A^, etc.) that must be set 
within the OSM algorithm. Some of these parameters, called secondary pa- 
rameters^ do not greatly affect the evolution, and can be set the same for all 
Nissl-stained images. The secondary parameters are as follows: is set to 
20, so that for a typical 80 fim circumference of a neuron cell body, neigh- 
boring control points are 3 /xm, or roughly 4 pixels away from each other. 
and S define the training that depend on the typical shapes of the object of 
interest, in our case a neuron. We build these parameters by creating a sam- 
ple of 100 ellipses, ranging linearly from an eccentricity of to 0.4, a simple 
representation of the average shape of neuron cell bodies. To speed up the 
evolution, we allow for different "time" steps and Mumford-Shah parameters 
in Eqs. (12][6). In Eq. 12, dt^dtc and A^Ac. In Eq. |6| dt^dtu and A^A^ 



In this schema, dtc^ dtu^ and A^^ can be set as secondary parameters which 
do not need to change for any of the pictures. We set dtc = 100, dtu = 0.05, 
and Xu = 1- 

In addition to the secondary parameters, there are two primary parame- 
ters which greatly affect segmentation, and must be determined empirically: 
the energy ratio a between Ems and £'c, and the energy parameter Ac within 
the Ems term. 

Because the active contour algorithm described above was designed for 
generic object recognition, the algorithm itself (in addition to the free pa- 
rameters) can be "tuned" for the task of finding dark elliptical features that 
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are overlapping or relatively close to each other on a lighter background. We 
adjust the above algorithm in a simple way to accommodate overlapping: if 
fij — Uij{t) > near and inside the given control point, the contour is "leak- 
ing" out to find the edge of another feature next to it. We therefore multiply 
this control point's contribution to Ems by a free parameter r] greater than 
1. Here, 77 is a secondary parameter, and is set to 1.5 for all images. 
We now discuss each step in ANRA. 



2.3 Step I: Image Acquisition 

We test ANRA on Nissl-stained tissue samples of seven young adult (6.4-11.8 
years; mean 8.5 years) and seven aged (24.7-32.9 years; mean 30.1 years) fe- 
male Rhesus monkey subjects that were part of an ongoing study of the 



effects of aging on cognitive function (Cruz et a/., 2004). For each subject 



eight (4 from each of 2 sections) gray scale (1-256) 512x512 pixel images with 
1.5 pixels/micron resolution (^150 neurons per image) were taken from area 
46, layer 3 of the prefrontal cortex in the ventral bank of sulcus principalis. 
3 subjects had appreciable differences in image quality between the two sec- 
tions, therefore the total number of different subject /image-qualities is 17. 
Fig. [2|3 shows 12 of the 17 subject /image-qualities. 



2.4 Step II: Segmentation Training 



All images are normalized as described in Sec. |2.1[ Out of each of the 17 
subject /image-qualities, one image is randomly selected as a training image. 
The digital image is marked for neuron cell bodies by an expert observer 
who "paints" sets of pixels over the neurons using a small graphical pro- 
gram. Different objects can share pixels, or overlap, but the sets exist as 
separate entities even if there is an overlap. We designate these sets of pixels 
created by an expert observer the training segments. The training segments 
will be compared to computer segments from the OSM output. The manual 
identification is relatively quick (2-4 seconds per neuron), and does not re- 
quire a model image, ie: no feature overlap ( Lin et al| |2005D . Furthermore, 



the cell marking method creates knowledge of the extents of each cell body 
as viewed by an expert observer, independent of and unbiased to our segmen- 
tation procedure. This information is saved and used repeatedly for multiple 
training runs as needed, and does not have to be repeated for the same image 
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if different training parameters are checked (Lin et a/., 2005; de Solorzano 



et a/., 1999). 



We next determine the values of the primary parameters a and the 
two primary free parameters which greatly affect the segmentation. We find 
that there is significant loss in functionality when a is outside the [10~^, 10~^] 
range and Ac is outside the [1,5] range. We therefore search this space of 
a and Ac by comparing the resulting computer segments to the training 
segments. A training segment is "found" if the computer segment shares more 
than 70% of the pixels with the training segment (Fig. [9]). The set called the 
final OSM parameters, denoted a* and Ac*, is the set that correctly identifies 
95% or more training segments. The (a*, Ac*) values are then recorded and 
used for the rest of ANRA. 

The OSM with the correct primary parameters (a*, Ac*) identifies 95% or 
more of neurons in the images, but it also identifies other non-neuron objects, 
such as staining errors, glial cells, and improper coverings of neurons. To 
separate neurons from non-neurons, computer training is performed. 

First, we compare the (a*, Ac*)-parameter OSM computer segments to 
the training segments. Each computer segment is either placed in the neu- 
ron segment category or non-neuron segment category based on whether the 
segment mutually overlaps any training segment (Fig. [9]). Second, each seg- 
ment is represented by seven segment properties v = ('^i, '^2, '^7). The seven 
segment properties were chosen to be the most salient measures of identifying 
neurons within an image. For the calculations of the segment properties, we 
denote the total number of pixels within the segment as Ac and the total 
number of pixels within the contour as \C\. The properties are based on the 
optical density of the original image fij as well as the square of the magni- 
tude of the image gradient \^fijf- The segment properties are presented in 
Table I. is a sum over all of the pixels within the segment area, is a 
sum over the edge pixels of the segment circumference, Tc is the location of 
the center of the segment, and r^j is the location of the pixel (ij). 



Using the WEKA machine learning toolkit ( Witten fc Frank[ 2005), we 
assess the following machine learning algorithm's ability to discriminate be- 
tween neuron property vectors |v ^, vj", . . .| and non-neuron property vectors 



{^1 ^ ^2 ^ ••• }' the 1-rule classifier (Holte, 1993), naive Bayes classifier (John 



fc Langley] |1995 ), support vector machine classifier (Piatt, 1998), nearest 



neighbor classifier (Aha & Kibler, 19911), decision tree classifiers (Quinlan 



|l993|), Bayes net and multi-layer perceptron ( [Witten fc Frank , 2005). The 
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cost between Type 1 errors (marking a non-neuron property vector as a neu- 
ron) and Type 2 errors (marking a neuron property vector as a non-neuron) 
is scanned by tuning the cost ratio term in the training algorithm. A strati- 
fied cross-vahdation evaluation for various cost ratios (3:1, 2:1,. ..,1:3) creates 
a receiver operator characteristic (ROC) curve ( |Duda et a/. , 2001) for each 
training method (Fig. [To|). The Multilayer Perceptron (MLP) using a single, 
4-node hidden layer, has the best ROC curve, as it provides the highest per- 
centage of neuron property vectors identified and the smallest percentage of 
non-neuron property vectors incorrectly identified. MLP is therefore chosen 
as the main training method for ANRA. 



2.5 Step III: Application 

Automatic neuron recognition is now applied on an unlimited number of 
other images that are normalized and similar in morphology to the training 
images. The steps are as follows: 

1. The OSM with the primary parameters (a*, Ac*) is performed on the 
new image. 

2. The properties v are calculated for each computer segment. 

3. A cost ratio is selected by the user. 

4. All computer segments deemed non-neurons by the MLP are discarded. 

5. For any two remaining computer segments that mutually overlap by 
more than 70%, the computer segments with the smaller probability of 
being a neuron (as determined by the MLP) is discarded. 

The (x,y) centers, sizes, and shapes of the remaining computer segments 
are the final result of ANRA. 



2.6 Comparison method 



A semi-automatic method (semi-auto) was used in prior neuron density maps 
correlation studies (ICruz et al\ 120051). In the semi- auto method a combina- 



tion of computer software and human intervention for each image is employed 
to identify neurons. Because the amount of human intervention scales with 
the number of images analyzed, the semi-auto method represents a standard 
with which we evaluate our completely automated recognition method. 
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2.7 Density Map Method and Microcolumnar Strength 



We give a description of the density map method, as it is the main anal- 
ysis to be apphed to the results of ANRA. The density map method was 
initially described by Buldyrev et al (2000) and a more detailed description 
and validation was given by Cruz et al ( 2005 ) . The density map is a 2D rep- 
resentation of the density correlation function g(x,y), which uses as input the 
(x,y) locations of all neurons in the region of interest (ROI). This function 
g{x^y) is mapped to a two-dimensional gray scale image (density map) in 
which different shades of gray are proportional to the average local neuronal 
density. Thus, the density map quantifies the average neuronal neighborhood 
surrounding a typical neuron within the ROI. 

Operationally, the density map is calculated by first assigning indices 
{i = 1,2,3...A^) to all the neurons in the sample. Next, we center a grid of 
bins of size D over each neuron and count how many other neurons fall in 
each bin constructing one matrix of accumulated neurons m(x, y). We define 
g{x^y) = m{x^y)/N'D'2^ in which g{x^y) has units of an average density of 
objects at position (x^y). As an example, the density map would be uniform 
if locations of objects (neurons) are uncorrelated, but will show patterns 
when there are regular spatial arrangements between the objects. 

For the case of neurons forming microcolumns, their density map exhibits 
one central vertical ridge, sometimes accompanied by two less pronounced 
parallel neighboring ridges. For this study, we are interested in the micro- 
column strength aS, which is extracted from the density map by taking the 
ratio of the neuronal density within the average microcolumn to the average 
neuronal density (Cruz et a/., 2005). For the same images, S is calculated 



using ANRA (x, y) locations as well as semi-automatic (x, y) locations, and 
the results are compared. 



3 Results 

For each of 17 subject /image-qualities, an evaluation image is randomly se- 
lected from the remaining images and marked for neuron cell bodies by the 
expert. The evaluation image is used as a "gold standard" to assess the 
accuracy of ANRA and the comparison methods. A total of 2448 "gold stan- 
dard" neurons are analyzed, for an average of 144 neurons per subject /image- 
quality. For each of the two recognition methods (semi-auto and ANRA), we 
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compare the method's identified neurons to the "gold standard", and retrieve 



the foUowing numbers (Fig. 11): 



a = number of correctly identified neurons , (13) 

b = number of non-neurons incorrectly identified as neurons , (14) 

and 

c = number of non-identified neurons . (15) 

To compare methods for the different subject /image-qualities, we define the 
following normalized metrics: 

A = . 100 , (16) 

and 

B = — ^ . 100 . (17) 

a + c ^ ^ 

A is the percent of correctly identified neurons ( "true positives" ). B is the 

percentage of non-neurons that are incorrectly identified as neurons ("false 
positives" ) . 



The results are shown in Fig. [T2j The semi-auto method is characterized 
by one {A^ B) set. Because of the ability to adapt the cost ratio as described 
in Sec. 2.5| ANRA is shown at 7 different ratios (3:1, 2:1, 1:1, 1:2, 1:3, 1:5, and 



1:10), ranging from very selective, to no selectivity, creating an "adapted" 
ROC curve. Since each point is an average of the 17 subject /image-qualities, 
the error bars show the standard deviation of the spread for both A and B. 
We choose the 1:2 cost ratio for further analysis because it is at the infiection 
point of the "adapted" ROC curve, and it has the closest average (A,B) to 



that of semi-auto. Table II and Fig. 13 a shows the individual results for 



each subject /image-quality for the semi- auto method and the ANRA with 
1:2 cost ratio. Fig. [T3| 3 shows an example of semi- auto and ANRA points 
compared to the gold standard. 

The results show that ANRA has a significantly higher A value of recog- 
nition (P-value: 0.002) and a similar B value of recognition compared to the 
semi-auto method. 



We also compare microcolumnar strength S (Sec. 2.7) using the (x^y) 
locations from both ANRA and semi-auto methods of neuron identification 
for the entire image database of rhesus monkey subjects as described in 



Sec. 2.3, 14,000 neuron locations were used, for an average of 1000 neuron 
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locations for each subject. We find significant correlations between micro- 
columnar strength measurements of the ANRA and semi-auto methods of 
neuron recognition (Fig. 14). This shows that ANRA has the ability to find 
significant changes in advanced neuron spatial arrangements within different 
subjects, and can therefore be applied to large datasets where manual or 
semi-auto recognition are not viable. 



4 Discussion 

In the present work we introduce a method called an Automated Neuron 
Recognition Algorithm (ANRA) which uses a combination of image segmen- 
tation and machine learning to retrieve neuron locations within digitized im- 
ages of Nissl-stained Rhesus monkey brain tissue. Despite challenges, such as 
overlapping of neuron cell bodies and the presence of glial cells and artifacts in 
the tissue, we demonstrate that ANRA has a significantly better recognition 



capability than a semi- auto method (Cruz et a/., 2005) which requires expert 



manual intervention for each image. ANRA's recognition quality is combined 
with computational efficiency, resulting in recognition of '^lOO neurons per 
minute using a standard personal computer. Consequently, large numbers 
of neuron locations can be retrieved, spanning considerably larger brain re- 
gions than ever before. Furthermore, because ANRA is capable of efficiently 
extracting neuron locations from durable and commonly used Nissl-stained 
tissue, it can potentially be applied to vast stores of archival material existing 
in laboratories and research collections around the world. 

Such a large dataset of (x,y) neuron locations will allow for a variety 
of systematic analyses that have previously not been possible. The ability 
to identify every neuron in entire sections of the brain will allow for both 
global and local analyses of neuron numbers, glial cell numbers, regional cell 
densities, and local variations in cell densities. Also, as was shown in the Re- 
sults section, studies of microcolumnarity or other spatial features of cortex, 
including spatial inter-relationships among neurons and glia using autocorre- 
lation and cross-correlation, are possible. Lastly, ANRA also allows for less 
obvious applications, including the investigation of the spatial network of the 
brain using the neuron locations as nodes. None of these studies are possible 
with the elegant sampling methods of modern stereology. 

We highlight the need for large datasets of neuron locations (10^ — 10^) 
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in comparative studies proposed in the Introduction and defined in Sec. |2.7 



Generally, the goal of a comparative study is to find a statistically signifi- 
cant diflFerence in a measured quantity (i.e., microcolumn strength) due to 
a change in an independent variable (age, species, sex, disease state, etc.). 
In the case of a ID correlation between nearest neighbors or the 2D mi- 
crocolumnar analysis, the neuron locations are used to create ID and 2D 
histograms, respectively. The number of neurons must be high enough to 
resolve the effect of the independent variable above random noise of the his- 
togram. Buldyrev et al ( 2000D showed that for a resolution of interest (seeing 
3% changes between 10 micron bins), ^ 10^ neuron locations are needed in 
the comparative study of microcolumnarity. For the same resolution in a 
ID correlation comparative study, such as nearest-neighbor distances, only 
^1000 neurons are needed (Schmitz et a/., 2002). For a given bin size, the 
theoretical calculation shows that the required number of neurons scales as a 
power of dimensions that are being correlated. Thus, automatic recognition 
becomes critical in higher dimension correlations. As an example we con- 
sider a 30 subject study of neuron spatial arrangement using 10^ neuron 
locations, making 100 different measurements of 1000 neurons each through 
a certain layer across several Brodmann regions. The semi-automatic ap- 
proach, which allows for acquisition of 10 neurons per second, would take 83 
human hours to complete. Comparatively, ANRA could complete the same 
task in 24 hours on 20 Intel P4 processors with less than 1 hour of preparation 
time. 

ANRA has a further advantage of reducing experimental drift. Specifi- 
cally, in terms of human bias, the "criteria" for neuronal identification will 
necessarily differ between different observers that are often required for a huge 
analysis extending over months to years, while ANRA's criteria, once estab- 
lished from the training algorithm, remains constant. Furthermore, ANRA's 
criteria will not be subject to the kind of experimental drift that can occur 
over time when one observer manually identifies thousands of neurons over 
a period of weeks to months. 

Recently, there have been advances in level set methods to recognize over- 
lapped cell nuclei (Lin et a/., 2007; de Solorzano et a/., 1999). The recognition 
challenges with Nissl-stained tissue are far greater than the challenges using 
confocal microscopy using fiuorescence. Lin et al ( 2007[ ) show how neurons 
and glia cells completely separate into two regions of parameter space using 
only two parameters (texture and intensity) of the identified segmentations. 
If plotted in a similar way, no two parameters that we consider (size, inten- 



16 



sity, texture, gyration, edge vs. area, etc.) would yield such a separation. 
Thus, in a Nissl-stained tissue visualized by optical microscopy, the param- 
eterized method of Cremers et al (2000), which, by design, overcomes the 
challenges of noisy images and missing boundary data (Sec. 2.2.2), is most 
efficient. 

Our results suggest that the ANRA method is performing as maximal 
efficiency: when a second expert's marks are compared with the gold standard 
on the same Nissl-stained image, the performance (^4 = 88 ± 5%) is not 
significantly higher than ANRA's performance (yl = 86 ± 5%). 

Although there are 10 free parameters within the algorithm, only two 
of them called the primary parameters must be explored to find the correct 
values for proper segmentation. These primary parameters are automatically 
found in the OSM parameter search during training. The other eight free 
parameters, which we call the secondary parameters^ can be fixed for the 
general task of identifying elliptical features within noisy images with missing 
boundary data, thereby solidifying them for the broadly applicable problem 
of neuron recognition in all Nissl-stained tissue. For a given morphological 
feature of interest, once a small set of representative images have been trained 
to, the training and parameters can be reused, due to the normalization of 
images of different quality. This setup will allow for the study of large areas 
of montaged images, or large datasets of hundreds of slides, all with the same 
training. Furthermore, the free parameters and training can be adapted for 
identification of other types of neurons, glial cells, etc. 

Lastly, because of the modular nature of the method (Fig. [T]), it will be 
relatively easy to replace partial aspects of the overall algorithm by consid- 
ering advances in recent published work. For example, [Tscherepanow et~ar 



(2006) independently developed a method to identify living cells that uses a 
larger set of training properties that is reduced with principle/independent 
component analysis, and Costa & BoUt (2006) has applied advanced pattern 
matching to the identification of neuron cell bodies in Nissl-stained tissue. 
By replacing the respective aspects of ANRA with such methods, the ideal 
overall identification algorithm can be found for not only the recognition 
of neuron cell bodies, but also the recognition of other objects of scientific 
interest, for example living cells or glial cells. 
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5 Software 



The ANRA software is available at http : / /physics . bu . edu/~ainglis/ 
flNRA/. 
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Table I 





description 


equation 


Vl 


segment area 


Ac 


V2 


average optical density (/) 




V3 


variance of optical density 


A^ Yl {fij ~ f) 


V4 


radius of gyration of optical density 


A^ ~ '^clfij 


V5 


segment edge length ( C ) vs. segment area 


\C\/A, 


Vq 


average gradient of segment edge 




vr 


average change in gradient of segment edge 
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Table II 





semi-auto 




ANRA 




# 


A(%) 


B(%) 


A(%) 


B(%) 


1 


81 


13 


82 


11 


2 


71 


14 


84 


7 


3 


82 


14 


91 


21 


4 


79 


15 


78 


4 


5 


90 


43 


92 


30 


6 


65 


3 


85 


16 


7 


76 


18 


87 


21 


8 


83 


15 


88 6 


9 


76 


12 


93 


15 


10 


79 


10 


85 


23 


11 


73 


6 


77 


7 


12 


82 


4 


88 11 


13 


92 


44 


84 


6 


14 


90 


26 


95 


16 


15 


80 


20 


80 


11 


16 


77 


23 


91 


28 


17 


75 


7 


86 


17 


avg. 


80±7 


17±12 


86±5* 


15±8 
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Figure 1: A schematic diagram showing processes involved in the Automated 
Neuron Recognition Algorithm (ANRA). The schematic describes the two 
main steps of the algorithm: training and application. Rectangles denote 
parameters that pass through the algorithm. Ovals, such as the OSM, are the 
computational parts of the algorithm, which can have images, segmentations, 
and parameters as their inputs and outputs. 
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(b) 




Figure 2: Challenges of automated neuron recognition, (a) 20x micrograph 
(scale bar: 50/xm) of a typical section showing the difficulties of separating 
neurons from glial cells and other artifacts in Nissl-stained tissue: 1. capil- 
laries, and unidentified material, 2. large glia (astrocytes), 3. glial as light 
as neurons in some cases, 4. neurons overlapped by glia (oligodendrocytes), 
5. neurons overlapped by other neurons, 6. multiple neurons and glial over- 
lapped, (b) lOx micrograph examples showing varying image quality. The 
highlighted micrograph is selected as an "ideal" contrast to be used in image 
normalization. 
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Figure 3: (a) Preprocessing "normalizes" the images so that they every image 
has the same background and foreground average optical densities, thereby 
removing the challenge of varying image type within Nissl-stained tissue. 
This is done by mapping optical density values of non-ideal images to an 
ideal image so that the average foreground and background averages are the 
same. The graph shows the optical density ranges of the ideal and non-ideal 
images (0..255), and a Bezier curve that passes through 4 points: (0,0), the 
background and foreground averages of the ideal and non-ideal images, and 
(255,255). (b) Examples of image normalization. 
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Figure 4: Steps of the overall segmentation method (OSM). (a) Over-marking 
the image with a hexagonal grid of points that lay on the thresholded fore- 
ground and center points of a traditional watershed segmentation. Points 
within 5 pixels are combined to avoid redundancy, (b) Active contour seg- 
mentation: using each starting location found in (a), a segmentation (clus- 
tering) process is performed within a small region of the image to find one 
possible neuron cell body. This process is then repeated for each starting lo- 
cation until all starting locations are exhausted, (c) The final set of computer 
segments^ shown in different solid colors, is the output of the OSM. 
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control point n 




Figure 5: (a) The contour C is described by the control points ri, r2, r^v. 
(a) Quadratic Bezier curve B{t) is defined for control point n using the control 
points r^_i, r^, and r^+i. The points r/^ are halfway between r^_i and r^. 
The equation for the contour is B{t) = (1 — t)^r/^ + 2t(l — t)r^ + r/^+it^, t = 
0..1. The equations guarantee that at the points r/ the curve is continuous 
and smooth. Combining several Quadratic Bezier curves creates a quadratic 
B-spline contour. An example with 5 control points is presented in (b) which 
shows how the B-spline contour moves when one control point (ri) moves, 
(c) Contour C (white pixels) with 20 control points (single black pixels) that 
is overlaying the image. 
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Figure 6: Schematic drawing showing the relative energies of Ems and Ec 
for the same image (shown as gray) and four different contour shapes (shown 
as black loops). The first three cases are examples of improperly fit contours 
with a high overall energy E = Ems + olEc. The last case is an example of 
an optimal contour minimizing the overall energy. 
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Figure 7: Examples of the active contour movement within the image during 
the OSM segmentation phase. The number of control points (black dots) 
is 20. The B-spline contour is white. The contour starts at a location de- 
termined by the over-marking step of the OSM. (a) Evolution of the energy 
terms Ems and Ec (b) Contour evolution after 0, 4, 8, 30, 48, 60, and 120 
steps. When a local minima is reached, the contour no longer moves, and 
the points internal to the contour are saved. 
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Figure 8: Control point movement based on Ems follows Eq. [12} The terms 
e+ and e~ are Ems (Eq. [s]) integrated over the single line of pixels right 
outside (e+) and right inside (e~) of the contour centered around each control 
point n. Tlx and are the x and y components of the outer normal vector 
of C at the control point. 




Figure 9: Two segments represent the same object when they mutually share 
more than 70% of their pixels. The two segments in (a) do not pass the 
required criteria because neither segment overlaps the other by more than 
70%. The two segments in (b) do not pass the required criteria because 
only one segment overlaps the other by more than 70%. Only in (c) does 
the required overlap occur. This analysis is used when computer segments 
are compared to "gold standard" training segments and either designated 
a neruon or non-neuron, and during the overlap deletion phase, when the 
segment with the highest probability of being a neuron is selected among all 
overlapped segments. 
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False Positive % 

Figure 10: Receiver operating characteristic (ROC) curve for each training 
method evaluated. It is seen that the Multilayer Perceptron (MLP) has the 
best ROC curve - the highest percentage of neuron property vectors identi- 
fied with the smallest percentage of non-neuron property vectors incorrectly 
identified. MLP is chosen as the main training method for ANRA. 
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Figure 11: Venn diagram showing the relative quantities for evaluating the 
quality of a neural recognition method. The bold black line separates neuron 
from non-neuron objects in the image. The dotted area shows the objects 
that are identified by a method. The method correctly identifies most of the 
neurons (a), but misses some neurons (c) and identifies some non-neurons 
as neruons (b). Using the quantities a,b, and c, standardized percentages of 
neuron vs. non-neurons can be calculated. 
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Figure 12: Results of ANRA. The semi-auto method is characterized by one 
{A^ B) set. Becasue of the abihty to adapt the cost ratio as described in 
Sec. |2.5| ANRA is shown at 7 different ratios (3:1, 2:1, 1:1, 1:2, 1:3, 1:5, and 
1:10), creating an "adapted" ROC curve. Since each point is an average of 
the 17 subject /image-types, the error bars show the standard deviation of 
the spread for both A and B. 
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Figure 13: (a) Individual results for 17 subject /image-types for the semi- 
auto method and the ANRA (with 1:2 cost ratio), (b) Recognition results 
for the semi-auto method (left) and the ANRA method (right) for example 
subject /image-quality #1 (Table II). Dark green: gold standard marks that 
match with the method. Blue: gold standard marks that DO NOT match 
with the method. Light Green: method points that match with gold standard 
points. Pink: method points that do not match with gold standard points. 
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Microcolumnar Strength 
(Microcolumn Density / Global Density) 




Figure 14: Comparison of microcolumnar strength measurement (S) using 
the {x^y) locations from both ANRA (with 1:2 cost ratio) and semi-auto 
methods of neruon identification. A total of 14,000 neuron locations were 
used, for an average of 1000 neuron locations for each subject (plot point). 
Both the neuron density and microcolumnar strength show significant corre- 
lations of ANRA with the semi-auto method. 
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